Equilibrating temperature-like variables in jammed granular subsystems 
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Although jammed granular systems are athermal, several thermodynamic-like descriptions have been pro- 
posed which make quantitative predictions about the distribution of volume and stress within a system and 
provide a corresponding temperature-like variable. We perform experiments with an apparatus designed to gen- 
erate a large number of independent, jammed, two-dimensional configurations. Each configuration consists of 
a single layer of photoelastic disks supported by a gentle layer of air. New configurations are generated by 
alternately dilating and re-compacting the system through a series of boundary displacements. Within each con- 
figuration, a bath of particles surrounds a smaller subsystem of particles with a different inter-particle friction 
coefficient than the bath. The use of photoelastic particles permits us to find all particle positions as well as 
the vector forces at each inter-particle contact. By comparing the temperature-like quantities in both systems, 
we find compactivity (conjugate to the volume) does not equilibrate between the systems, while the angoricity 
(conjugate to the stress) does. Both independent components of the angoricity are linearly dependent on the 
hydrostatic pressure, in agreement with predictions of the stress ensemble. 

PACS numbers: 45.70.-n, 64.30.-t, 83.80. -Fg 



Granular materials are a collection of discrete, athermal 
particles. In the absence of an external driving force, these 
materials relax into a mechanically stable jammed state and 
cannot move into another configuration since thermal fluctu- 
ations are negligible (TJ. While these materials are therefore 
inherently non-equilibrium, preparing a configuration with a 
strict protocol nonetheless yields different microscopic states 
with the same, reproducible volume 0. Edwards proposed 
that the system volume (a conserved quantity) could be used 
to write a granular density of states, a corresponding entropy, 
and a temperature-like variable conjugate to the volume J3). 
However, a complete granular statistical mechanics should de- 
scribe the distribution of contact forces as well as the volumes. 
Subsequent theoretical advances have proposed that a stress- 
based ensemble B WT1 is likely required for a full treatment. 

In the Edwards ensemble, the volume V plays a role anal- 
ogous to that of energy in equilibrium statistical mechanics. 
A granular temperature, dubbed the compactivity, is defined 
as X = (dS/dV)^ 1 , and has been successfully measured in 
models |[T2l [13) , simulations |[T4l [T31 , and experiments |[T6l - 
|2D . Similarly, the stress ensemble considers force and torque 
constraints on individual particles, and writes the density of 
states as a function of the stress-tensor S = ?iifij> where 
the fij are the vectors pointing from the center of each parti- 
cle to its contacts, and /y is the corresponding contact force. 
The conjugate variable is then a tensorial temperature known 
as the angoricity, and is defined to be A = (dS/dJ^)^ 1 . 

A minimal test of such temperature-like variables is to con- 
sider whether they obey the zeroth law of thermodynamics. In 
experiments and simulations, the compactivity ||20ll has pre- 
viously been shown to be equal in different parts of the same 
packing, and in different packings generated with the same 
particles under identical conditions. Simulations show this is 
also satisfied by the angoricity QH). However, no test has 
been made of whether two dissimilar systems can equilibrate 
either X or A. We provide such a test in a real granular sys- 



tem subject to isotropic compression, and find that while the 
compactivity fails this simple test, the angoricity equilibrates 
in a temperature-like way. 

Our experiments are conducted on a bi-disperse granular 
monolayer of photoelastic disks resting on a nearly friction- 
less surface provided by a thin layer of pressurized air. The 
assembly of particles is comprised of an inner subsystem and a 
larger bath which differ only in the inter-particle friction coef- 
ficient (see Fig.[T|. Starting from a dilute state, the monolayer 
is bi-axially compressed by outer walls in a series of short 
steps. At some global volume fraction <P, the system jams and 
for all further steps the pressure on the system increases. Fi- 
nally, the walls re-dilate to permit large scale rearrangements 
before the next series begins. By repeating this protocol many 
times, we generate an ensemble of configurations for which 
we record particle positions and calculate contact forces using 
methods similar to 1221 1231 . With this information, we calcu- 
late the compactivity and angoricity for both the bath and the 
inner subsystem. 

In the canonical volume ensemble [3|, the probability of 
finding a system with volume V and compactivity X is pro- 
posed to be given by a Boltzmann-like distribution 



V(V) = W e -v/x 



(1) 



where the density of states is fl(V) is defined for an ensemble 
of jammed configurations, and the partition function is Z(X). 
The stress ensemble similarly proposes a Boltzmann-like dis- 
tribution 



7>(E) = ^ e-^ ( g /^ 
Z(A) 



(2) 



for the stress-moment tensor S; the angoricity A is therefore 
also a tensor. 

To calculate either X or A, we use two methods: the 
method of overlapping histograms ||7] |20l |24) and the 
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fluctuation-dissipation theorem (FDT) liTfjl [TTl |2D . The ra- 
tio of V(V) between two systems is exponential in V and is 
given by 



-Pi(V) 
V 2 (V) 



Z(X 2 ) 
Z(X t ) 



(3) 



By taking the logarithm of this ratio, one obtains a term lin- 
ear in V, where the coefficient is the difference in the inverse 
temperatures. This method determines 1/X up to an additive 
constant: l/X — > 1 jX + Cx ■ The FDT method also provides 
a relative measurement. Using the measured variance (8V 2 ) 
of V(V), we compute 
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to obtain values of X, also up to a constant. The calculation 
of A utilizes equations analogous to Eqn. [3] and Eqn. [4J the 
tensorial aspects will be discussed in more detail below. Each 
of these methods is used separately on both the subsystem and 
the bath, in order to test for equilibration. 

Our experimental apparatus is shown to scale in Fig.[T] The 
granular monolayer consists of 1004 bi-disperse photoelastic 
(Vishay PhotoStress PSM-4) disks with a thickness « 3.1 mm 
and diameters d$ — 11.0 mm and d^ = 15.4 mm, in equal 
concentrations. The particles are supported on a thin layer 
of air provided by a steady flow of pressurized air through 
a porous polypropylene sheet with a nominal pore size of 
120 /im. This minimizes the effect of friction between the 
particles and the surface, but does not otherwise cause sig- 
nificant dynamics once the system is jammed. The sheet is 
leveled (particles do not drift to one side) and flat (particles do 
not cluster). The system consists of an outer bath Nb = 904 
and an inner subsystem Ns = 100. Particles in the bath have 
a friction coefficient [Ib ~ 0.8, while particles in the inner 
subsystem are wrapped with a thin layer of PTFE tape with a 
/is < 0.1. 

Images of the particle positions, photoelastic images for 
measuring vector contact forces, and identification of the sub- 
system particles are recorded with three separate images cap- 
tured by a single CCD camera located above the apparatus 
(see Fig. [TJ)). Particle positions are identified using a white 
light image (see Fig.[T]:), from which the centers are detected 
with an accuracy of « O.Olris using a Hough transform. The 
photoelastic images (see Fig. [T]l) are captured using reflective 
photoelasticity, in which the silvered back side of each parti- 
cle reflects polarized light back to the camera. Photoelasticity 
allows for the numerical determination of the normal and tan- 
gential forces at each contact point, as required to measure S. 
Similar to the methods pioneered by 11221 l23l . we minimize 
the error between the observed and fitted image of the parti- 
cle using a non-linear least-squares optimization. Details and 
source code are available for download at ll25ll . The third im- 
age is taken using black-light illumination to identify the sub- 
system particles, which are tagged with ultraviolet-sensitive 
ink (see Fig.[lJ). The subsystem comprises all low-/i particles 
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FIG. 1. Color online. Schematic of apparatus showing (a) two walls 
bi-axially compressing an array of disk-shaped particles composed 
of an outer subsystem (black, high /i) and an inner subsystem (red, 
low [£) and (b) reflective photoelasticity on air-floated particles. Light 
shines from green LEDs through a linear polarizer (P) a wavelength- 
matched quarter wave plate (Q) before entering the photoelastic ma- 
terial. A mirrored surface on the bottom of each particle reflects light 
back through the particle. A second quarter-wave plate and linear 
polarizer are mounted on the camera to resolve the photoelasticity. 
Three images of each configuration are recorded: (c) unpolarized 
white light for locating particle positions, (d) polarized green light 
showing isochromatic fringes for calculating contact forces and (e) 
an ultraviolet light for identifying the low-/i particles. 



which are Voronoi neighbors with at least one other particle in 
the subsystem. 

The particles are confined within a square region (maxi- 
mally 50 x 50 cm) imposed by two stationary walls positioned 
by stepper motors, as shown in Fig.[lJ. The system is initially 
in a dilute, un-jammed state, with the global volume fraction 
$ < 0.6. The two walls bi-axially compress the system by a 
series of small steps of constant size (A$ — 0.0009, equiva- 
lently Ax = 0.3 mm or 0.02 r^). With each step of the wall, 
the three images are recorded, and data is collected over a 
series of volumes corresponding to 0.775 < $ < 0.805, giv- 
ing 30 different volumes for each compression cycle. Steps 
continue until the gradient squared of the force image ll26ll 
indicates a pressure threshold has been reached; this reduces 
the risk of particles buckling out of plane. The walls then re- 
dilate to the dilute state, and the particles are then rearranged 
while maintaining subsystem continuity; this protocol is re- 
peated 100 times. 

During the compression phase of each quasi-static cycle, 
we observe the percolation of force chains at a value $ perc . 
As the system is further compressed beyond this point, the 
contact forces grow in strength and the average number of 
contacts per particle increases. For the set of 100 cycles, this 
threshold occurs over a range 0.782 < ^ per c < 0.792, where 
the width of the distribution is indicative of finite size effects 
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FIG. 2. Color online, (a) Volume histograms, V(V), for $ = 0.776 
(A), 0.784 (■), and 0.802 (•) with m = 48. (b) A semi-logarithmic 
plot of the ratio each histogram with respect to the $ = 0.784 dis- 
tribution, i.e. Vi(V)/Vi=2(V). (c) The inverse compactivity given 
by Eqn.[3]plotted as a function of the inverse volume fraction where 
\ib are shown as black • and (is are red ♦. Large/small symbols 
denote jammed/un-jammed configurations, respectively. Errorbars 
shown are uncertainties in V(V) and propagated through the calcu- 
lation. The inverse compactivity given by the FDT method (Eqn.|4|, 
is shown with the solid line for comparison, (d) The ratio of number 
of jammed/un-jammed configurations recorded at each $. 



Il27l l28l . The ratio of un-jammed to jammed systems at a 
given <£> is shown in Fig. [2]l We define random loose packing 
as $flxp = ($ P erc) ~ 0.787 as the center of this distribution. 

We calculate the distribution of local volumes V(V m ) over 
clusters of size to, using the sum of individual radical Voronoi 
volumes obtained from the Voro++ software |29|. Each clus- 
ter is defined as the to — 1 nearest neighbors surrounding a 
central particle. For m = 1, V(V m ) has two distinct peaks 
which correspond to small and large particles ifTSl . With in- 
creasing cluster size, the bimodal aspect of V(V m ) disappears, 
but even for large cluster sizes (to > 100), the distribution 
remains asymmetric and non-Gaussian [18 |. In Fig. [2^, we 
show V(V m ) for three values of $ with to = 48; the value of 
to is large enough so that V(V m ) does not show any features 
arising from bi-dispersity. 

In Fig.|2j5, we show the ratio Vi(V)/Vj(V) where the ref- 
erence system j is <E> = 0.784. In practice this can be done 
with any two systems so long as there is sufficient overlap 
between their histograms. As the ratio of ~Pi(V)/Vj(V) is 
well-approximated by an exponential in V, the compactiv- 
ity can be calculated using Eqn. [3] The inverse compactiv- 
ity, 1/X, is also calculated using FDT using Eqn. |4] where 
the integrand is approximated using a third order polynomial. 
Each method determines 1 jX only up to an additive constant, 
which is adjusted so that Xrlp = oo. In Fig. |2j;, the in- 
verse compactivity is shown for both the bath (1 /Xb) and the 
subsystem (1/Xg). We find good agreement between X(§) 



given by the overlapping histogram method and by the fluctu- 
ation dissipation theorem. In addition, for 4 < to < 50, we 
observe X to be approximately independent of to. However, 
we find that the compactivity of the bath is not equal to that of 
the subsystem ^ Xs($)), even considering adjust- 

ments of the additive constant. This represents a failure of the 
zeroth law for X. 

We can take further advantage of the accessibility of both 
jammed and un-jammed states within in the center of the 
range of explored $. While the Edwards ensemble is not de- 
fined for un-jammed systems, we can nonetheless carry out 
the histogram analysis as performed on the jammed systems. 
In this regime, we find that the ^(Vm) histograms cannot dis- 
tinguish between the jammed and un-jammed states. Further- 
more, the measured values of X decrease continuously from 
above Qrlp to below; this is an undesirable characteristic. 

The stress ensemble also provides a Boltzmann-like distri- 
bution in the stress, as given in Eqn. [2] In the case of friction- 
less grains, the angoricity A is a scalar due the off-diagonal 
components in £ being zero. In any real granular system, 
friction is present and a shear-free state is not readily obtained. 
Therefore, S is a symmetric tensor with non-zero off-diagonal 
components and can be reduced to two independent compo- 
nents related to the pressure and shear stress. The pressure 
angoricity A p and the shear angoricity A T are conjugate to 
(T p = (<ti +(i2)/2 and the a T — (a\ — 02V2, respectively lf8l . 
where a± t 2 are the principal stresses. The average hydrostatic 
pressure per particle in the system is given by T = Tr T,/N . 
Both A p and A T are obtained using the method of overlapping 
histograms (analogous to Eqn. [3]) and the FDT (analogous to 
Eqn. [4]). With each method, A is calculated up to an addi- 
tive constant so that 1/A —> 1/A+ Ca, where Ca satisfies 
A 00 as T — >• 00. 

In Fig. [3^, the local distribution of pressure V(a p ) is shown 
for to = 8 on configurations over a range 0.0006 < T < 
0.0025 Nm. The ratio Vi(a p ) /Vj(a p ) is exponential in <j p 
(see Fig. ^p, similar results for a T not shown), as required 
by the stress ensemble analogue of Eqn. [3] In addition, we 
observe that the variance of a p is proportional to to, which 
is consistent with 5 being an extensive entropy (see Fig.|3j;). 
We are therefore able to measure both the pressure angoricity 
A p the shear angoricity A T using their corresponding distri- 
butions, shown in Fig.|3j; as a function of T. We find that A PtT 
are independent of to for to > 3, as also observed in sim- 
ulations 10 E), and that values obtained from the histogram 
method (points) and the FDT method (solid line) are in ap- 
proximate agreement. Finally, we find that for either the shear 
or compressional angoricity, the values measured in the bath 
and in the subsystem are equivalent, signifying the angoricity 
is equilibrating between the subsystems. 

Nonetheless, the values of A T and A p do not match each 
other, with the shear angoricity growing faster as a func- 
tion of r. We find the angoricity is given by A = b T for 
both pressure angoricity and shear angoricity, where b p = 
0.153±0.004and6 r = 0.450±0.020, respectively. For a two- 
dimensional frictionless shear-free system, the stress ensem- 
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FIG. 3. Color online, (a) Distribution of a v where m — 8 and T = 
0.0007 (T), 0.0010 (■), 0.0015 (•), and 0.0024 Nm (A). A semi- 
logarithmic plot of the (b) ratio V% (a p ) /Vj (o>) where the reference 
system j is T = 0.0015 Nm. The pressure angoricity Ap and shear 
angoricity A T are shown as a function of V where the results using 
overlapping histograms for [ip and \jls are shown as black o and are 
red <0>> respectively. The solid line is the angoricity calculated using 
FDT The gray dashed lines provide a visual reference of the slopes 
0.15 and 0.45, respectively. Inset: The scaled variance (8<j^) of the 
Vj (<7 P ) distribution, as a function of the cluster size m. 



ble predicts b p — 0.5 at the isostatic point QUO- Above the 
isostatic point, the stress ensemble predicts b p to be a function 
of the average contact number. The disagreement between the 
frictional and frictionless values of b p implies friction signifi- 
cantly affects the density of states. 

We have measured both compactivity X (conjugate to vol- 
ume in the Edwards ensemble), and angoricity A (conjugate to 
the stress tensor in the stress ensemble), in a laboratory gran- 
ular system using particle-scale characterizations. While we 
found that while the value of X calculated using the overlap- 
ping histogram method was consistent with the value found 
using the fluctuation-dissipation theorem, it failed to equili- 
brate between non-identical systems, making it a poor state 
variable. A similar failure is likely behind previous measure- 
ments by Schroter et al. 11711 . in which two granular materials 
with different frictional properties, prepared using the same 
protocol, were found to have different globally-measured val- 
ues of X. In contrast, we observed that the temperature-like 
variable A does successfully equilibrate between a subsystem 
and bath with dissimilar inter-particle friction coefficients, as 
would be required in order to have a valid zeroth law. More- 
over, we find agreement with the prediction that angoricity 
should scale linearly the hydrostatic pressure [8 |. These suc- 
cesses make angoricity a promising state variable for frictional 
granular systems. 

One downside to using angoricity as a state variable, partic- 
ularly in experiments, is that its calculation requires the deter- 



mination of both normal and tangential forces. While there 
has been a long history of measuring normal forces at the 
boundaries of granular systems [ 30-34 1, particle-scale mea- 
surements have seen more limited development. Outside of 
photoelastic particles such as those used here, measurements 
typically exist only for normal forces, whether the systems are 
frictional (tangential forces are neglected) |35 , 36 1 or friction- 
less BTHSa. 



It is possible to understand the success of the stress ensem- 
ble over the Edwards (volume) ensemble by considering the 
underlying physics behind the conserved quantities in each. 
Under Newton's third law, forces and torques must be strictly 
balanced at each force contact, while volume is merely con- 
strained globally. As a result, our subsystem differed from the 
bath not only in the measured X, but more conventionally in 
the mean local volume fraction. 

In fact, the full canonical Edwards ensemble [6] unifies 
the volume and stress ensembles, where the density of states 
depends on both V and S, and it has recently been ar- 
gued BOl |4D that the two should not be considered sepa- 
rately. The classic phenomenon of Reynolds dilatancy ll42ll 
under which shear induces a bulk expansion similarly sug- 
gests that such a coupling is important. Nonetheless, we ob- 
served here that angoricity can be independently equilibrated, 
and future experiments should more fully investigate the rela- 
tionship between ensembles, the relative importance of shear 
and compression, and the role of friction on the density of 
states. 
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